Self-consistent electronic structure of a d x i_ y i and a d x 2_ y 2 + id xy vortex 
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We investigate quasiparticle states associated with an isolated vortex in a d-wave superconductor 
using a self-consistent Bogoliubov-de Gennes formalism. For a pure d x 2_ y 2 superconductor we find 
that there exist no bound states in the core; all the states are extended with continuous energy 
spectrum. This result is inconsistent with the existing experimental data on cuprates. We propose 
an explanation for this data in terms of a magnetic-field-induced d x 2_ y 2 +id xy state recently invoked 
in connection with the thermal conductivity measurements on Bi2Sr2CaCu2 0s. 



Behavior of superconductors in a magnetic field has 
been traditionally at the center of the condensed matter 
research because of the rich variety of fascinating phe- 
nomena in such systems and also because of its consider- 
able technological importance. With the advent of high- 
T c cuprates the subject became even more interesting 
owing to the realization that their order parameter is un- 
conventional, exhibiting most likely a d x 2_ y 2 symmetry. 
A number of novel effects associated with the existence 
of low-lying quasiparticle excitations near the gap nodes 
have been predicted . While there exists experimen- 
tal evidence for each of these effects (4f](|, it is by no 
means conclusive and an active debate continues. Sev- 
eral fundamental issues remain unresolved. Among the 
most interesting is the electronic structure of magnetic 
vortices, where a coherent physical picture is still lack- 
ing despite the number of theoretical and experimental 
investigations. Notably the nature of quasiparticle core 
states in d-wave superconductors remains poorly under- 
stood. These states are of considerable importance as 
they impact on various static and dynamic properties of 
the mixed state such as the structure of the flux lattice, 
pinning and flux-flow resistance. 

It was established more than three decades ago by 
Caroli, de Gennes and Matricon that discrete quasi- 
particle states exist localized in vortex cores of con- 
ventional s-wave superconductors. These states are la- 
beled by the angular momentum quantum number /z, 
and the low-lying eigenvalues are E y ~ ^(Aq/Ef) with 
/i = 1/2,3/2,..., Ao the bulk gap and Ep the Fermi 
energy. This prediction has been later confirmed in de- 
tail by numerical computations and by experiments 
on NbSe2 In a d x 2_ y 2 superconductor the situation 
becomes considerably more complex owing to the non- 
trivial structure of the gap function which vanishes along 
the four nodes on the Fermi surface. While the s-wave 
bound states can be intuitively understood by drawing 
an analogy to a simple quantum-mechanical problem of 
a particle in the cylindrical well of the radius £ ~ vf/ttAq 
and height Ao, a suitable analogy in a (i-wave case would 
involve a potential well whose radius and height depend 
on the polar angle 8, with A(0) vanishing along the four 
diagonal directions causing to diverge. Under such 
circumstances one would expect low-lying states to be ex- 



tended along the node directions, rather than localized, 
unless the mixing of the core levels into the continuum is 
prevented by some higher symmetry. The existing exper- 
imental work, however, appears consistent with localized 
core states with large energy spacing [p|jic| |l2[| . 

The problem of a vortex in a d x 2_ y 2 superconductor 
has been considered theoretically by a number of au- 
thors [p~3| — but the detailed nature of quasiparticle core 
states has not been addressed. Very recently Morita, 
Kohmoto and Maki studied this problem using an 
approximate version of the Bogoliubov-de Gennes (BdG) 
theory. As pointed out in a subsequent Comment [ p)| |, 
this approximation improperly neglects an essential ele- 
ment of the physics of d x 2_ y 2 superconductors and yields 
unphysical results. 

In this article we present, for the first time, a fully self- 
consistent numerical solution of the continuum BdG the- 
ory for a single isolated vortex in a d-wave superconduc- 
tor. In agreement with the above qualitative argument 
we find that for a pure d x 2_ y 2 gap there exist no truly 
localized core states pofl . We find low lying states that 
are strongly peaked in the core but have tails along the 
gap node directions which do not appear to decay to zero 
far from the core. The energy spectrum of these states 
becomes continuous in the limit of infinite system size. 
This is consistent with their extended nature, but incon- 
sistent with the experimental finding of a large gap, 
Eq « Ao/5, to the lowest core state in YBa2Cu307_«5. 
We propose that one can reconcile the theory with ex- 
periment by assuming a parity and time reversal sym- 
metry violating d x 2_ y 2 + id xy state which has been re- 
cently invoked to explain the thermal conductivity data 
on Bi2Sr2CaCu20g in finite magnetic fields M,E2|. Such 



a state is fully gapped and for sufficiently large d xy com- 
ponent it can support true bound states in the core of 
the nature similar to the s-wave case. 

The BdG equations for a d-wave vortex have been pre- 
viously solved numerically on the lattice Jl3|jl4|]. In this 
work we adopt a continuum formulation of the problem, 
which has been used to study conventional s-wave vor- 
tices B. Such formulation allows us to directly exploit 
various symmetries of the system and as a result we are 
able to consider much larger systems with a reasonable 
spectral resolution. The BdG equations for a d-wave su- 
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perconductor can be written as E3] 

H e u(R) + J drA(R - r/2, r)v(R - r) = Eu(R), 
-7t!>(R)+ /drA*(R-r/2,r)u(R-r) =J3«(R). (1) 



Here 7i e is the single electron Hamiltonian, which we 
assume to have a simple free particle form 



-A 

c 



(2) 



and A(R,r) = V(r)(c|(R — r/2)cj.(R + r/2)) is the or- 
der parameter which is a function of the center-of-mass 
and the relative coordinate R and r. In terms of the 
quasiparticle amplitudes [u n ,w n ] the order parameter is 
self-consistently determined from the gap equation 

A(R, r) = V(r)J2'K(R - r/2)<(R + r/2) 

+ u„(R + r/2)<(R - r/2)] tanh(£„/2T), (3) 

where n labels the eigenstates of (|l|) and the prime means 
the usual restriction to energies smaller than a cutoff 
scale fi c . V(r) is the pairing interaction which we as- 
sume to have the following model form: 



V(r) = V 5(r) + V l9 (<p)-6{r-a), 
a 



(4) 



with r = (r, ip) . The potential (^) consists of a contact re- 
pulsion term Vq > 0, meant to model the on-site Coulomb 
repulsion between holes, and a short range attractive part 
Vi < 0, necessary to establish superconducting pairing. 
In order to model the lattice structure of cuprates we 
allow this attraction to be angle dependent with 



g((f) = (1 - e) cos 2 (2(/j) + esin 2 (2^), 



(5) 



where e is a parameter used to tune the relative strength 
of this model's analog of nearest and second nearest 
neighbor attraction. The nature of superconducting in- 
stability in this model depends on the dimensionless pa- 
rameter rj = akp. In the absence of magnetic field the 
ground state has a <i-wave symmetry for 7/ > 1.9 and is 
extended s-wave for rj < 1.9. This is consistent with the 
phase diagram of the related lattice model, where d-wave 
is known to be stable only close to half filling pifl . Fur- 
thermore, in the d- wave regime, the ground state is a pure 
d x 2_ y 2 for e = 0, a pure d xy for e = 1, and a d x 2_ y 2 +id xy 
admixture for < e < 1. More generally one can also 
allow for anisotropy in the single particle Hamiltonian 
TL e to account for the band structure effects in cuprates. 
We find that such anisotropies do not modify the phase 
diagram significantly and have only secondary effect on 
the vortex core structure. In the following we therefore 
limit ourselves to the simple form (Q) and we furthermore 
neglect the vector potential A, as appropriate for the ex- 
treme type-II cuprates . The magnitude of the contact 



repulsion Vq has little effect on the physics in the d-wave 
regime For simplicity we therefore take Vq = 0. 

The form of the pairing interaction ([|) implies that 
the order parameter depends on the relative coordinate r 
only through its polar angle tp. It is convenient to expand 
it in terms of the 2D angular momentum eigenstates, 

A(R, r) = A(R, 9 W )=Y, e- ip9 e^A pl iR), (6) 
P .i 

where R = (R, 9) . In such a representation the inte- 
ger p characterizes winding of the superconducting phase 
around the vortex and I specifies the orbital state of the 
Cooper pair. Thus, for instance, the dominant order pa- 
rameter near a singly quantized d x i^ y i vortex consists 
of an equal superposition of p = 1, I = ±2 components: 
A(R, 6; if) = 2e~ i0 cos(2^)A 1)±2 rR). 

We solve the BdG equations (|l|) numerically on a disk 
of the radius Rq with a suitably chosen initial order pa- 
rameter and we then iterate Eqs. ([l]) and (|J) until self- 
consistency is achieved. Following Gygi and Schluter 
Q we expand the quasiparticle amplitudes in the basis 
spanned by the eigenfunctions of TL e : 



[«(R),w(R)] 



/i,m 



<$>» m {R)[u firm v fj 



(7) 



Here $ Mm (i?) = [V2/i?o^+i( a A"n)] ^Km^o), ^( z ) 
is the Bessel function of order /x and a Mm is the m-th 
zero of Jf_i{z). The integro-differential equation (|l|) thus 
becomes an eigenvalue problem with an infinite matrix 
which we truncate at large values of and m, and di- 
agonalize using a standard LAPACK subroutine. In the 
s-wave case this matrix is block diagonal in the angular 
momentum /i and the resulting radial problem can be 
solved for each fj, separately. The crucial new element in 
the d-wave case is that, due to the complicated structure 
of the pairing term (^) , the angular momentum channels 
remain coupled and one has to solve the full 2D problem. 
This essential feature was ignored in fll8|| . 

To model a pure d x i_ y i case we choose 77 = 2, e = 0, 
fl c /E F = 0.3, Vx/Ep = 1.6 and T = 0. We assume 
the initial order parameter of the form Ax,±2{R) = 
A tanh(i?/£ ), where £0 = VF/nAd is the coherence 
length p5j. The above parameters imply A^ = 0.26Ep 
and £0 = 2.5kp 1 . In the process of iterating Eqs. (|l|) 
and (^) various other components A p i(R) appear in the 
self consistent solution reflecting the spatial anisotropy of 
the dominant d x 2_ y 2 component and nucleation of vari- 
ous sub-dominant order parameters near the core. Fig. 
[l] shows the leading components of the self-consistent or- 
der parameter. The d x 2_ y 2 order parameter relaxes to its 
bulk value over a distance ~ 20k p 1 , which is much larger 
than £o- in fact, by varying the coupling strength V\ in a 
wide range, we find that for coherence lengths £o^f ^5 10 
the order parameter profile (rescaled to its bulk value) 
no longer depends on £0 and retains the universal shape 
showed in Fig. [l| We also find that this profile is not 
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FIG. 1. Dominant order parameter components for sys- 
tem size RokF = 120. In terms of conventional notation we 
have \d x 2_ 2] oc (Ai,_ 2 + Ai, 2 ), \d xy \ oc (Ai l _ 2 — Ai,a) and 



'A_i,o + e 



-3i<p 



very well fitted by the usual tanh(i?/£') function for any 
£' but is more consistent with algebraic ~ 1/r 2 relax- 
ation. Fig. [I] also shows sizable components with s-wave 
symmetry of the form predicted by the Ginzburg-Landau 
theory We have explicitly verified that, in agreement 
with |H , they will form four satellite vortices at a distance 
~ 26kp 1 form the origin. 

We now investigate the nature of quasiparticle core 
states by computing the generalized inverse participation 
ratios, defined as pal, 



a n = (R k F ) 2 



(KI 4 )s + (KI 4 )s 
((Kl 2 ) s + (kl 2 ) s ) 2 



(8) 



where (. . .) s — j . . . dR. As a function of increasing sys- 
tem size i?o this quantity approaches a finite constant for 
an extended state and grows as i? 2 , for a state localized 
within the characteristic length £t <C Rq. Fig. || shows a n 
as a function of energy for system sizes R^kp = 80, 120. 
Over the entire energy range the data for two differ- 
ent sizes behave in a similar way. If there existed lo- 
calized states in the core their corresponding a n would 
have increased by more than a factor of 2 between sizes 
RokF = 80 and 120. No such increase is observed. We 
carried out similar analysis for an s-wave case (same pa- 
rameters with rj = 1), where it is known that only local- 
ized states exist at low energies. Indeed, in this case a 
clear a n ~ Rq scaling was observed [E7| for E n < Aq. 
Fig. ^ further shows that, unlike in the s-wave case [||, 
there is no discernible pattern in the core energy leve s 
and their spacing decreases with increasing system size. 
Visual inspection of the amplitudes of these states (in- 
set Fig. ||) reveals that they are highly anisotropic with 
peaks near the core and tails running along the gap node 
directions which appear to saturate to a finite amplitude 
far from the core. Similar states are known to exist near 
a strongly scattering non-magnetic impurity [p8[ . 

Although not truly localized these core states should 
still give rise to enhanced tunneling conductance from 
the cores. In Fig. |^a we display this quantity computed 
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FIG. 2. Inverse participation ratios function of en- 

ergy. The arrow marks the maximum bulk gap Ad. Inset: 
amplitudes |u n (R)| and |u n (R)| of the lowest core state. 



from the present model. The closely spaced core states 
give rise to a broad peak centered near the zero energy. 
The peaks corresponding to individual states can still 
be resolved, but by studying the size dependence of such 
spectra we conclude that in the limit of infinite size these 
will form a continuum, similar to that reported by Wang 
and MacDonald in their lattice calculation |Q. While 
the fine details of the spectrum depend on parameters of 
the model as well as on the system size, its gross qual- 
itative features, and in particular the vanishingly small 
gap to the lowest state, remain remarkably robust for 
a wide range of parameters considered. These spectra 
are, unfortunately, in disagreement with the experimen- 
tal finding of a large ~ gap to the lowest core state 
0.0. It is in principle possible that the intragap spec- 
tral features found experimentally are in fact resonances 
rather that true bound states. Even so, we find no indica- 
tion of such resonant behavior in our calculation. Ref. Q 
also reports that the core states are spatially isotropic, in 
contradiction to strong anisotropies found by this work 
and by quasiclassical computations p^JbSl . 

It would thus appear that a model based on a pure 
d x 2_ y 2 order parameter is in all aspects inconsistent with 
the published experimental data. We therefore propose 
an explanation in terms of a mixed d x 2_ y 2 + id xy state, 
which may be formed in cuprates in finite fields at low 
temperatures |2l],|2^]. Such a state is fully gapped and 
we therefore expect the quasiparticle core states to re- 
cover some of their localized s-wave character. Within 
the present model we study such a scenario by taking 
e = 0.25 (with other parameters unchanged). This results 
in a d x 2_ y i + id xy state with d xy ~ 0.17d x 3_ y a, roughly 
consistent with Laughlin's prediction [^] for T = and 
B = 6T. The order parameter distribution near the core 
remains qualitatively similar to that reported in Fig. [l], 
but now with unequal magnitudes of Ai^ and Ai i _2- 
Carrying out the analysis of participation ratios we find 
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FIG. 3. Tunneling conductance g(-E, R) at the center of the 
vortex (solid line) and far from the vortex (dashed line), both 
spatially averaged over a circular ring 2kp 1 wide. Dashed and 
dash-dotted vertical lines mark the maximum and minimum 
bulk gap respectively; RokF = 120. 



that the core states below the minimum gap d xy are truly 
localized and nearly isotropic, in agreement with exper- 
iment. Fig. shows the corresponding tunneling con- 
ductance, which is now qualitatively consistent with ex- 
periment in that there exists a finite gap to the lowest 
core state which is independent of the system size. This 
gap is however still too small to account for the experi- 
mental result; we find that the agreement becomes better 
for larger d xy component. Note that within the present 
model the d xy component is brought about by change in 
the pairing interaction, rather than magnetic field and it 
would therefore be present even at zero field. Although 
we expect that the qualitative features of the core states 
are insensitive to the origin of d xy , a more satisfactory 
model would correctly describe the transition to the time 
reversal breaking state as a function of field. This would 
require explicitly including the vector potential term A 
in the kinetic energy (Q); the work on this is in progress. 

Our results firmly establish the absence of localized 
vortex core states in superconductors with a pure d x 2_ y 2 
gap. We speculate that a sizable admixture of a subdom- 
inant order parameter is needed to reconcile theory with 
experiment. A magnetic- field-induced d xy component 
appears to be the most acceptable choice at present. This 
is by no means free of problems. A sufficiently large d xy 
component could be difficult to justify since on general 
grounds one would expect d xy j d x i_ y i to scale as the ratio 
of the corresponding critical temperatures p2[ | . Further- 
more, a transition to the fully gapped d x i_ y i + id xy state 
would be observable in the specific heat or muon-spin- 
rotation measurements of the penetration depth, but no 
such effect has been reported fH^. In the present con- 
text of the vortex core states, however, our proposition 
can be easily tested by measuring the core spectra at 
lower fields (or higher temperatures) . If the present spec- 



tra are indeed characteristic of a d x 2_ y 2 + id xy state a 
dramatic change should be observable at the transition 
to the pure d x 2_ y 2 state. 

The authors are indebted to A. V. Balatsky, A. J. 
Bcrlinsky and A. J. Millis for helpful discussions. This 
research was supported by NSF grant DMR-9415549. 
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